In this lab we’ll introduce the add-on package tmap, which can be used to produce much nicer maps than the default mapping functions in R. We’ll also look briefly at leaflet, a package for making interactive maps. Before starting the lab, make sure both of these packages are installed. Other useful ones include mapsf highcharter, mapview and mapdeck.
Remember to use RStudio’s script editor (top left panel) to enter your commands and then copy-paste these to the console when you are ready. This will keep a copy of the commands that you can save and return to.
Before starting, remember to create a working directory
(e.g. module15). Next download the files for today’s lab
from the Canvas page and move them to this directory (unzip the
compressed files). You’ll need the following files:
Note for the shapefiles, you will need to download and decompress the zip file for each one to get the shapefile and the associated extra files. If you are not sure about this, please let me know.
Now start RStudio, and change your working directory from the
[Session] menu > [Set working directory] > [Choose directory…].
Remember to check that R can see the files by running the
list.files() command.
Once you have done this, let’s load all the files we will need:
library(sf)
library(terra)
sf_use_s2(FALSE)
tmap makes thematic maps. It works in a very similar
way to ggplot2, by building a series of layers and map
geometries and elements. In general, we start by using
tm_shape() to identify the spatial object to be used, and
add geometries are added, including filled polygons, borders and
symbols. Finally, we can add legends, scale bars, etc. We’ll look at
three examples of working with this package here, but there are a lot
more details and examples on the [TMap website][tmapid].
We’ll first make a map with a mixture of vector layers, including census tracts, light rail tracks and stations. Let’s start by loading these:
## Census tracts for Salt Lake County with population density
tracts <- st_read("./slc_tract/slc_tract_2015.shp", quiet = TRUE)
## Salt Lake light rail tracks
lightrail <- st_read("./LightRail_UTA/LightRail_UTA.shp", quiet = TRUE)
## Salt Lake light rail stations
stations <- st_read("./LightRailStations_UTA/LightRailStations_UTA.shp", quiet = TRUE)
If we check the CRS, you’ll see that tracts has a
different one, so will need to be reprojected:
st_crs(tracts)$epsg
## [1] 4269
st_crs(lightrail)$epsg
## [1] 26912
tracts <- st_transform(tracts, st_crs(lightrail))
Now let’s start mapping things out with tmap. First,
let’s make a simple map showing the polygon outlines using
tm_borders(). This is the standard tmap
method: first define the object that you want to plot with
tm_shape, then define what you want to plot.
library(tmap)
tm_shape(tracts) +
tm_borders()
The function tm_fill() will then fill these using one of
the variables in the tract data set. We’ll use the population density
(density). Note that this automatically adds a legend
within the frame of the figure:
tm_shape(tracts) +
tm_fill("density") +
tm_borders()
The color scale can be changed by setting the palette
argument in tm_fill(). This includes ColorBrewer and
viridis scales and there are a set of methods todefine the intervals in
the color palette. For example, to use the ‘Greens’ palette with
percentile breaks. We also change the color of the tract boundaries to
lightgray:
tm_shape(tracts) +
tm_fill("density", palette = "Greens", style = "quantile") +
tm_borders("lightgray")
You can reverse the color scale by prepending a
-. This
uses a reversed magma palette from the viridis set, with a
continuous palette:
tm_shape(tracts) +
tm_fill("density", palette = "-magma", style = "cont") +
tm_borders("lightgray")
If you want to see the full set of palettes that you can use, install the tmaptools package and run the following code (you may also need to install shinyjs):
tmaptools::palette_explorer()
Let’s now add another layer. We’ll overlay the light rail track on
this map. As this is a different object in R (lightrail),
we first need to use tm_shape to indicate that we are using
it, then add tm_lines() to show the routes. We’ll make the
lines a little thicker with lwd and change the type to
dashed.
tm_shape(tracts) +
tm_fill("density", palette = "Greens", style = "quantile") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2, lty = 'dashed', col = "darkorange")
Note that if you use a varible name for the color, it will thematically map the lines (I’ve removed the polygon fill to make this clear):
tm_shape(tracts) +
#tm_fill("density", palette = "Greens", style = "quantile") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 4, col = "ROUTE")
We can now add the stations (using
tm_shape() to select the
correct object to plot):
tm_shape(tracts) +
tm_fill("density", palette = "Greens", style = "quantile") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2) +
tm_shape(stations) +
tm_dots(size = 0.25, shape = 23)
Now let’s add some details to the map. We’ll add a graticule before
the polygons are plotted, a compass and a scalebar. We’ll also use the
tm_layout function to add a title and move the legend to a
new position. (This function has a lot of options for changing your map
- it’s well worth checking the help page.)
tm_shape(tracts) +
tm_graticules(col = "lightgray") +
tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2) +
tm_shape(stations) +
tm_dots(size = 0.25, shape = 23) +
tm_compass(position = c("left", "bottom")) +
tm_scale_bar(position = c("right", "top")) +
tm_layout(main.title = "Salt Lake County Light Rail",
legend.outside = TRUE,
legend.outside.position = c("left"))
The map can be clipped to smaller regions by using the
bbox argument in the first tm_shape()
function. In this code, we’ll extract a set of census tracts that
correspond to downtown Salt Lake. We’ll then use the the bounding box of
these to crop the map (I’ve removed the fill). We’ll also add the
station name using tm_text():
tracts_sub = tracts %>%
dplyr::filter(TRACTCE %in% c(102500, 102600, 114000))
tm_shape(tracts, bbox = st_bbox(tracts_sub)) +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2) +
tm_shape(stations) +
tm_dots(size = 0.25, shape = 23) +
tm_text("STATIONNAM", ymod = -1, bg.color = "white", size = 0.8)
Finally, tmap allows you to make interactive maps
fairly simply with the tmap_mode() function. Setting this
to view will make an interactive map, to plot
will make a static map (like the ones we have made so far). Here we’ll
remake the full map as an interactive map:
## Set interactive
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(tracts) +
tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile", id = "TRACTCE") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2) +
tm_shape(stations) +
tm_dots(size = 0.25, shape = 23)
## Symbol shapes other than circles or icons are not supported in view mode.
We’ll reset back to static maps for the rest of this exercise:
tmap_mode("plot")
## tmap mode set to plotting
In this example, we’ll map out the climate data from western North
America. We’ll go a bit faster in this example and not explain every
step. First we’ll need to convert the wna_climate data
frame to an sf object:
wna_climate <- read.csv("./WNAclimate.csv")
wna_climate <- st_as_sf(wna_climate,
coords = c("LONDD", "LATDD"),
crs = 4326)
Make a quick thematic map of average july temperature:
tm_shape(wna_climate) +
tm_symbols(col = "Jul_Tmp")
Next, we’ll download country outlines and river centerlines from
[Natural Earth][natEarthID] using the rnaturalearth
package (you will need to install this). This will download layers from
the Natural Earth website - there are a variety of these available,
including country and state boundaries, physical features (rivers,
lakes, etc) and cultural features (place names, roads, etc). There are
three levels of resolution: 1:10m, 1:50m and 1:110m. To download a
layer, you need to specify the resolution (scale) and the
name of the layer. The full set of available layers can be found
[here][rneid]
library(rnaturalearth)
countries50 <- ne_download(scale = 50, type = "admin_0_countries")
rivers50 <- ne_download(scale = 50, type = "rivers_lake_centerlines", category = "physical")
Now we can put this all together as a series of layers. We’ll make the color palette continuous and use a red-blue palette.
tm_shape(countries50, bbox = st_bbox(wna_climate)) +
tm_borders() +
tm_shape(rivers50) +
tm_lines("lightblue", lwd = 2) +
tm_shape(wna_climate) +
tm_symbols(col = "Jul_Tmp", palette = "-RdBu", alpha = 0.75,
style = "cont", title.col = "degC")
## Warning: The shape rivers50 contains empty units.
Finally, we’ll make two maps (july temperature and annual
precipitation), and use tmap_arrange() to make a two panel
plot:
m1 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
tm_borders() +
tm_shape(rivers50) +
tm_lines("lightblue", lwd = 2) +
tm_shape(wna_climate) +
tm_symbols(col = "Jul_Tmp", palette = "-RdBu", alpha = 0.75,
style = "cont", title.col = "degC")
m2 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
tm_borders() +
tm_shape(rivers50) +
tm_lines("lightblue", lwd = 2) +
tm_shape(wna_climate) +
tm_symbols(col = "annp", palette = "BuPu", alpha = 0.75,
style = "cont", title.col = "mm/yr")
tmap_arrange(m1, m2)
## Warning: The shape rivers50 contains empty units.
## Warning: The shape rivers50 contains empty units.
## Warning: The shape rivers50 contains empty units.
## Warning: The shape rivers50 contains empty units.
In the final example, we’ll use the Landsat data with tmap.
## Landsat images for California
filenames <- paste0('./rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- rast(filenames)
## Places for California
ca_places <- st_read("./ca_places/ca_places.shp", quiet = TRUE)
The main function to map raster data is tm_raster().
We’ll remake the NDVI layer that we made in lab 12 and plot this.
b4 <- landsat[[4]]
b5 <- landsat[[5]]
ndvi <- (b5 - b4) / (b5 + b4)
plot(ndvi, col=rev(terrain.colors(10)), main = "NDVI")
tm_shape(ndvi) +
tm_raster()
We’ll update this a little by clamping the NDVI values at a lower limit of 0, changing the color palette and making it continuous. We’ll also add the California places added as an overlay.
tm_shape(clamp(ndvi, lower = 0)) +
tm_raster(palette = "Greens", style = "cont", title = "NDVI") +
tm_shape(ca_places) +
tm_borders(lwd = 2) +
tm_layout(main.title = "Landsat 8 (2017/06/14)",
legend.position = c("left", "top"),
legend.bg.color = "white", legend.bg.alpha = 0.7)
A second function (tm_rgb()) allows the creation of
color composites. We’ll use the red, green and blue bands to make the
true color composite. The maximum reflectance value in these files is
1.0, so we need to set this for scaling.
tm_shape(landsat) +
tm_rgb(r = 4, g = 3, b = 2, max.value = 1)
## Warning: Raster values found that are outside the range [0, 1]
When we made this figure earlier, it was a lot brighter. This is
because the previous function ‘stretched’ the range of values for each
layer. We can mimic that here by first using the stretch()
function from terra, then repeating the plot:
landsat_stretch = stretch(landsat, minv = 0, maxv = 1,
minq = 0.01, maxq = 0.99)
tm_shape(landsat_stretch) +
tm_rgb(r = 4, g = 3, b = 2, max.value = 1)
leafletThe tmap mode function allows you to quickly make an
interacitve map. The resulting map is based on the well-known Leaflet
javascript library, which offers a much greater range of flexibility in
making these maps. Fortunately, there is an R library
(leaflet) that allows access to the Leaflet API
allowing you to build quite complex maps from R. We’ll remake the map of
Salt Lake County and the light rail routes using this now. First install
the library (install.packages("leaflet")), and load it (and
a helper library for working with html formats):
library(leaflet)
library(htmltools)
leaflet operates in a similar way to
ggplot2 and tmap in building a map by
layers. Layers are added using the pipe operator that we encountered in
an earlier lab (%>%). There are three basic ingredients
that are required:
leaflet() functionHere’s a really simple example, where we use the default tiles from OpenStreetMap and drop a marker on the University of Utah:
m <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m # Print the map
You can also try changing the Provider - the source of the background map. For example, this will use ESRI map (the map is not shown here but will appear in RStudio when you run this):
m <- leaflet() %>%
addProviderTiles("Esri.WorldStreetMap") %>%
addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m
ESRI imagery:
m <- leaflet() %>%
addProviderTiles("Esri.WorldImagery") %>% # Add default OpenStreetMap map tiles
addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m
You can find a full set of these providers here: https://leaflet-extras.github.io/leaflet-providers/preview/
We’ll now add the census tracts and light rail sf
objects. First, these need to be reprojected to WGS84 long-lat (EPSG
code 4326) to work with Leaflet.
tracts_ll <- st_transform(tracts, crs = 4326)
stations_ll <- st_transform(stations, crs = 4326)
lightrail_ll <- st_transform(lightrail, crs = 4326)
Now we can add these to a dark basemap. First the tracts with
addPolygons, then the routes (addPolylines)
and stations (addMarkers). Note that the stations include a
label taken from the station name column in the
sf object. This wil appear when you hover your cursor over
a marker.
leaflet() %>%
# add a dark basemap
addProviderTiles("CartoDB.DarkMatter") %>%
# add the polygons of the clusters
addPolygons(
data = tracts_ll,
color = "#E2E2E2",
opacity = 1, # set the opacity of the outline
weight = 1, # set the stroke width in pixels
fillOpacity = 0.2 # set the fill opacity
) %>%
addPolylines(
data = lightrail_ll
) %>%
addMarkers(
data = stations_ll,
label = ~htmlEscape(STATIONNAM)
)
Finally, we’ll add some extra details. First, we’ll set colors for each station by the light rail line it lies on. This is a little complicated - below we:
sf objectunique(stations_ll$LINENAME)
## [1] "University" "N/S TRAX" "Med Center"
## [4] "Intermodal Hub" "Airport" "Draper"
## [7] "Mid-Jordan" "West Valley" "Sugar House Streetcar"
line_pal <- RColorBrewer::brewer.pal(9, "Set1")
line_no <- as.numeric(as.factor(stations_ll$LINENAME))
stations_ll$LINECOL <- line_pal[line_no]
Next we’ll make a popup window for each station. This will appear
when the marker is clicked (as opposed to the label which appears when a
cursor hovers over the marker). This has to be formatted using html.
Here, we use R’s paste() function to stick together html
tags with the station name, address, etc. The <b>
tags make the labels bold, and the <br> tags add a
line break. You can do a lot more here, for example, including images or
URLs in each popup.
station_popup = paste0(
"<b>Station: </b>",
stations_ll$STATIONNAM,
"<br>",
"<b>Line Name: </b>",
stations_ll$LINENAME,
"<br>",
"<b>Park n Ride: </b>",
stations_ll$PARKNRIDE,
"<br>",
"<b>Address: </b>",
stations_ll$ADDRESS
)
With all this in place, we can rebuild our map. The main changes here are that we use a circle marker as this has an option to change color, and include the popup information in the marker.
leaflet() %>%
# add a dark basemap
addProviderTiles("Esri.WorldStreetMap") %>%
# add the polygons of the clusters
addPolygons(
data = tracts_ll,
color = "#E2E2E2",
# set the opacity of the outline
opacity = 1,
# set the stroke width in pixels
weight = 1,
# set the fill opacity
fillOpacity = 0.2
) %>%
addPolylines(
data = lightrail_ll
) %>%
addCircleMarkers(
data = stations_ll,
color = ~LINECOL,
label = ~htmlEscape(STATIONNAM),
popup = station_popup
)
There is no additional exercise for today, but you need to submit your R script or (preferably) a Markdown document with the code examples from the lab